!
!  MDark     Mean of a set of dark frames.
!
!  Copyright © 1998-2011 F.Hroch (hroch@physics.muni.cz)
!
!  This file is part of Munipack.
!
!  Munipack is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation, either version 3 of the License, or
!  (at your option) any later version.
!  
!  Munipack is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!  
!  You should have received a copy of the GNU General Public License
!  along with Munipack.  If not, see <http://www.gnu.org/licenses/>.

program mdark

!
! Computes mean of darks (masterdark) by the arithmetical mean 
! of all images (non-robust but fast) or by the robust mean (remove 
! defects on single frame but slow and memory consuming).
!

  use LibList
  use input
  use robustmean
  use stat

  implicit none

  ! Version:
  Character(len=*),parameter :: version = 'Mdark Version 0.3, (C) 1998-2010 F.Hroch (hroch@physics.muni.cz)'

  ! Initial parameters.
  !    Output image name:
  Character(len=32768) :: namedark='mdark.fits'

  !    Minimum of values
  Real :: minvalue = 0.0

  ! No. of image dimensions
  integer, parameter :: DIM = 2
  
  !-------

  Character(len=32768) :: sb,key,value,line,record
  Character(len=80) :: buf,dateobs,dateobs1
  character(len=8) :: TEMPERAT='TEMPERAT',DATE_OBS='DATE-OBS',EXPTIME='EXPTIME'

  Logical :: simple,extend,anyf!,exist
  Integer :: i,j,n, nobs
  Integer :: istat = 0, bpix, naxis, naxes(DIM)
  integer :: bitpix = 0
  Integer :: blocksize,pcount,gcount
  logical :: robust = .false.

! image dimensions
  Integer ::  nx = -1, ny = -1 
! image
  Real, dimension(:,:), allocatable :: DARK, SIG
  real, dimension(:,:), pointer :: ccd
  real, dimension(:), allocatable :: fbuf
  Real :: s, t, temp, stemp, e, etime,setime,maxvalue
!  real :: factor = 1.0
  type(Imagetype), pointer :: list,curr

  call InitList(list)
  curr => list

  do
     read(*,'(a)',iostat=i,end=20) record
     if( i /= 0 ) stop 'Failed to read of an input.'
  
     call input_rcomment(record,line)
!     write(*,*) trim(line),"<-",trim(record)

     select case(input_classify(line))

     case(input_key)

        call input_keys(line,key,value)
        i = 0
!     call input_read(key,value,i)
 !    if( i == opt_end .or. i == opt_err ) exit
 !    write(*,*) trim(key),':',trim(value)

        if( key == 'Output filename' ) then
           read(value,*,iostat=i) namedark
        endif

        if( key == 'Robust mean' ) then
           read(value,*,iostat=i) robust
        endif
             
        if( key == 'Output BITPIX' ) then
           read(value,*,iostat=i) bitpix
        endif

!        if( key == 'Product by factor' ) then
!           read(value,*,iostat=i) factor
!        endif
     
        if( key == 'Key TEMPERAT' ) then
           read(value,*,iostat=i) TEMPERAT
        endif
     
        if( key == 'Key DATE-OBS' ) then
           read(value,*,iostat=i) DATE_OBS
        endif
        
        if( key == 'Key EXPTIME' ) then
           read(value,*,iostat=i) EXPTIME
        endif

        if( key == '' .or. value == '' .or. i /= 0 ) then
           write(*,*) 'This record: ',trim(key),trim(value),' parse failed.'
        endif

!        write(*,*) trim(key)," ",trim(value)

!        write(*,*) robust
!        write(*,*) trim(namedark)

     case (input_data)

     ! input filename
!     read(*,*,end=20,iostat=i) sb
!     if( i /= 0 ) stop 'Bad format of input filename.'
!     write(*,*) trim(sb),':'

        read(line,*,iostat=i) sb
        if( i /= 0 ) then
           write(*,*) 'Failed to read of filename:',trim(line)
           stop 'Bad format of the input filename.'
        end if
     write(*,*) trim(sb),':'

     istat = 0
     ! input initialisation of the current file
     Call ftopen(25,sb,0,blocksize,istat)
     Call ftghpr(25,DIM,simple,bpix,naxis,naxes,pcount,gcount,extend,istat)
     if( istat /= 0 ) goto 666

     if( getno(curr) > 0 .and. (nx /= naxes(1) .or. ny /= naxes(2)) ) then
        write(*,*) "Dimensions of the image does not match. Skipping."
        call FTCLOS(25,istat)
        goto 665
     endif

     t = 0.0
     CALL ftgkye(25,TEMPERAT,t,buf,istat)
     if( istat /= 0 ) then
        write(*,*) TEMPERAT,' keyword not found.'
        istat = 0
     endif
     e = 0.0
     CALL ftgkye(25,EXPTIME,e,buf,istat)
     if( istat /= 0 ) then
        write(*,*) EXPTIME,' keyword not found.'
        istat = 0
     endif
     dateobs = ''
     CALL ftgkys(25,DATE_OBS,dateobs,buf,istat)
     if( istat /= 0 ) then
        write(*,*) DATE_OBS,' keyword not found.'
        istat = 0
     endif

     ! image - data input
     if( robust .or. getno(curr) == 0 ) then
        allocate(ccd(naxes(1),naxes(2)),stat=i) 
        if( i /= 0 ) then
           write(*,*) 'Not enought memory. Skipping.'
           call FTCLOS(25,istat)
           goto 665
        endif
     end if

     ! first image => initialisation
     if( getno(curr) == 0 )then
        nx = naxes(1)
        ny = naxes(2)
        
        Allocate(DARK(nx,ny),SIG(nx,ny),stat=i)
        if( i /= 0 ) stop 'Not enought memory for dark.'
        DARK = 0.0 
        SIG = 0.0
        dateobs1 = dateobs
     endif

     call FTG2DE(25,1,minvalue,naxes(1),naxes(1),naxes(2),ccd,anyf,istat)
     call FTCLOS(25,istat)
     if( istat /= 0 ) goto 666

     if( robust ) then
        Call AddItem (curr,array=ccd,dateobs=dateobs,exptime=e,temp=t)
     else

        ! compute sumation
        DARK = DARK + ccd
        SIG  = SIG  + ccd**2

        ! add image
        Call AddItem (curr,dateobs=dateobs,exptime=e,temp=t)
     endif

665  continue

  end select

enddo

  ! input of top line with init values
!  read(*,*,end=666,iostat=j) sb,i,e,bpix
!  if( j /= 0 )stop 'Bad format of input data.'
!  if( sb /= '' ) namedark = sb
!  if( i >= 0 ) robust = i == 1
!  if( e > 0.0 ) factor = e
!  if( bpix /= 0 ) bitpix = bpix

!  call InitList(list)
!  curr => list
!  do
!
!  enddo

20 continue
  nobs = getno(curr)
  if( bitpix == 0 ) bitpix = bpix
  
  if(  nobs == 0 ) stop 'No input image(s).'

  if( nobs > 1 ) then
     
     Allocate(fbuf(nobs),stat=i)
     if( i /= 0 ) Stop 'Not enought memory for a buffer.'

     ! computation of mean
     if( robust ) then

        Do j = 1,ny
           Do i = 1,nx
              
              curr => list
              do 
                 curr => GetNext(curr)
                 if( .not. associated(curr) ) exit
                 call GetItem(curr,array=ccd)
                 n = getno(curr)
                 fbuf(n) = ccd(i,j)
              enddo
              call rmean(fbuf(1:nobs),DARK(i,j), SIG(i,j))
              
           enddo
        enddo

     else

        DARK = DARK / nobs
        SIG = sqrt((max(SIG - nobs*DARK**2,0.0)/(nobs - 1.0)))
        
     end if

     ! computation of mean exposure time
     curr => list
     do 
        curr => GetNext(curr)
        if( .not. associated(curr) ) exit
        n = getno(curr)
        call GetItem(curr,exptime=fbuf(n))
     enddo
     if( robust ) then
        call rmean(fbuf(1:nobs),etime,setime)
     else
        call mean(fbuf(1:nobs),etime,setime)
     end if

     ! computation of mean temperature
     curr => list
     do 
        curr => GetNext(curr)
        if( .not. associated(curr) ) exit
        n = getno(curr)
        call GetItem(curr,temp=fbuf(n))
     enddo
     if( robust ) then
        call rmean(fbuf(1:nobs),temp,stemp)
     else
        call mean(fbuf(1:nobs),temp,stemp)
     end if
     
     deallocate(fbuf)

     ! computation of the mean of the standard deviation
     s = sum(SIG,mask=SIG>0) / count(SIG>0)

  else 

     ! nobs == 1
     DARK = ccd
     temp = t
     stemp = 0.0
     etime = e
     setime = 0.0
     s = 0.0

  end if

  ! product by factor
!  DARK = factor*DARK
!  s = factor*s
  ! the determination of s is perhaps only estimate in case of robust
  ! mean due to non-linear method

  ! cut-off range 
  if( bitpix > 0 )then
     maxvalue = 2.0**bitpix - 1.0
     where( DARK < minvalue )
        DARK = minvalue
     endwhere
     where( DARK > maxvalue )
        DARK = maxvalue
     endwhere
  end if

  ! remove output image, if exist
!  istat = 0
!  Inquire(file=namedark,exist=exist)
!  if( exist ) then
!     open(26,file=namedark)
!     close(26,status='DELETE')
!  endif

  ! save output image
!  istat = 0
  Call ftinit(26,namedark,1,istat)
  Call ftphps(26,bitpix,naxis,naxes,istat)
  Call ftpkys(26,'OBJECT','Mean Dark','',istat)
  call ftpkye(26,EXPTIME,etime,10,'Average of exposure times (seconds)',istat)
  call ftpkye(26,TEMPERAT,temp,10,'Average of camera temperatures (C)',istat)
  call ftpkys(26,DATE_OBS,dateobs1,'UTC of the first on input',istat)
  if( bitpix > 0 ) then
     call ftpkyj(26,'BSCALE',1,'',istat)
     call ftpkyj(26,'BZERO',2**(bitpix-1),'',istat)
  endif
  if( robust ) then
     write(buf,'(A,F0.1,A,F0.1)',iostat=i) &
          'Robust mean temperature = ',temp,'+-',stemp
  else
     write(buf,'(A,F0.1,A,F0.1)',iostat=i) &
          'Mean temperature = ',temp,'+-',stemp
  endif
  if( i == 0 ) call ftpcom(26,buf,istat)
  if( robust ) then
     write(buf,'(A,F0.1,A,F0.1)',iostat=i) &
          'Robust mean exposure time = ',etime,'+-',setime
  else
     write(buf,'(A,F0.1,A,F0.1)',iostat=i) &
          'Mean exposure time = ',etime,'+-',setime
  endif
  if( i == 0 ) call ftpcom(26,buf,istat)
  write(buf,'(a,f0.3)',iostat=i) 'Mean of sigmas = ',s
  if( i == 0 ) call ftpcom(26,buf,istat)

  if( nobs > 0 ) then
     write(buf,'(a,i0,a)',iostat=i) &
          'The image is a result of average of ',nobs,'  exposure(s).'
     if( i == 0 ) call ftpcom(26,buf,istat)
     call ftpcom(26,'The time of start and the temperature for each image used:',istat)
     curr => list
     do 
        curr => GetNext(curr)
        if( .not. associated(curr) ) exit
        call GetItem(curr,dateobs=dateobs,temp=t)
        write(buf,'(a, 1x, F0.1)',iostat=i) trim(dateobs),t
        if( i == 0 ) call ftpcom(26,buf,istat)
     enddo
  endif
  if( robust ) then
     call FTPCOM(26,'The robust mean has been used.',istat)
  else
     call FTPCOM(26,'The arithmetical mean has been used.',istat)
  end if
  call FTPCOM(26,version,istat)
!  write(*,*) istat

  if( istat /= 0 ) goto 666

  ! data 
  Call ftp2de(26,1,nx,nx,ny,DARK,istat)

  CALL ftclos(26,istat)
  if( istat /= 0 ) goto 666

  ! print some information on console
  write(*,*) 'Output image: ',trim(namedark)
  write(*,*) 'Dimension: ',nx,'x',ny,'    Mean of sigmas=',s

  if( .not. robust ) then
     deallocate(ccd)
  end if
  call DestroyList(list)
  deallocate(DARK,SIG)

666 continue 

  if( istat /= 0 )then
     write(*,*) 'FITSIO: following error(s) has been occured:'
     do
        Call ftgmsg(buf)
        if( buf == ' ' ) exit
        write(*,*) trim(buf)
     enddo
  endif

end program mdark
